function [new_weather] = SubSample(my_weather,time,alt, min_sigma, new_weather)
%*=+--+=#=+--      UAS Severe Weather Simulation Softeware        --+=#=+--+=#*%
%          Copyright (C) 2013 Regents of the University of Colorado.           %
%                             All Rights Reserved.                             %
%                                                                              %
%    This program is free software: you can redistribute it and/or modify      %
%    it under the terms of the GNU General Public License version 2 as         %
%    published by the Free Software Foundation.                                %
%                                                                              %
%    This program is distributed in the hope that it will be useful,           %
%    but WITHOUT ANY WARRANTY; without even the implied warranty of            %
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the             %
%    GNU General Public License for more details.                              %
%                                                                              %
%    You should have received a copy of the GNU General Public License         %
%    along with this program.  If not, see <http://www.gnu.org/licenses/>.     %
%                                                                              %
%           Jack Elston                                                        %
%           elstonj@colorado.edu                                               %
%                                                                              %
%*=+--+=#=+--                 --+=#=+--+=#=+--                    --+=#=+--+=#*%

		P = [];
		for x = my_weather.x
			for y = my_weather.y
				point = my_weather(time,x,y,alt);
				P = [P; point.u, point.v, point.w];
			end
		end

		C = P'*P;
		[V,D] = eig(C);

		sigma = 1-D(3,3)/sum(diag(D)) %FIXME, this seems incorrect
		%sigma = D(1,1)/sum(diag(D)) %FIXME, this seems incorrect
		p_bar = floor([mean(my_weather.x), mean(my_weather.y)]);

		if (sigma > min_sigma) || (size(P,1) > 1000)

%			[y,i] =	max(V(:,3));
%	    if i == 2
%				bounds.x0 = my_weather.x0;
%				bounds.xf = p_bar(1);
%				bounds.y0 = my_weather.y0;
%				bounds.yf = my_weather.yf;
%			else
%				if i == 1
%					bounds.x0 = my_weather.x0;
%					bounds.xf = my_weather.xf;
%					bounds.y0 = my_weather.y0;
%					bounds.yf = p_bar(2);
%				else
				  if size(my_weather.x,2) > size(my_weather.y,2)
						bounds.x0 = my_weather.x0;
						bounds.xf = p_bar(1);
						bounds.y0 = my_weather.y0;
						bounds.yf = my_weather.yf;
					else
						bounds.x0 = my_weather.x0;
						bounds.xf = my_weather.xf;
						bounds.y0 = my_weather.y0;
						bounds.yf = p_bar(2);
					end
%				end
%			end
		  bounds.z0 = alt;
		  bounds.zf = alt;

		  split_weather = Weather(my_weather.settings,my_weather.simdata,bounds);
			%plot(split_weather,split_weather.t0,split_weather.z0,'wind',1);
			%axis([0 19 0 19]);
			%axis([0 104500 0 104500]);
			%pause(0.1);
		  new_weather = SubSample(split_weather,time,alt,min_sigma,new_weather);

%	    if i == 2
%				bounds.x0 = p_bar(1)+1;
%				bounds.xf = my_weather.xf;
%				bounds.y0 = my_weather.y0;
%				bounds.yf = my_weather.yf;
%			else
%				if i == 1
%					bounds.x0 = my_weather.x0;
%					bounds.xf = my_weather.xf;
%					bounds.y0 = p_bar(2)+1;
%					bounds.yf = my_weather.yf;
%				else
				  if size(my_weather.x,2) > size(my_weather.y,2)
						bounds.x0 = p_bar(1)+1;
						bounds.xf = my_weather.xf;
						bounds.y0 = my_weather.y0;
						bounds.yf = my_weather.yf;
				  else
						bounds.x0 = my_weather.x0;
						bounds.xf = my_weather.xf;
						bounds.y0 = p_bar(2)+1;
						bounds.yf = my_weather.yf;
					end
%				end
%			end
		  bounds.z0 = alt;
		  bounds.zf = alt;

		  split_weather = Weather(my_weather.settings,my_weather.simdata,bounds);
			%plot(split_weather,split_weather.t0,split_weather.z0,'wind',1);
			%axis([0 19 0 19]);
			%axis([0 104500 0 104500]);
			%pause(0.1);
		  new_weather = SubSample(split_weather,time,alt,min_sigma,new_weather);

%		  bounds.x0 = p_bar(1)+1;
%		  bounds.xf = my_weather.xf;
%		  bounds.y0 = my_weather.y0;
%		  bounds.yf = p_bar(2);
%		  bounds.z0 = alt;
%		  bounds.zf = alt;
%
%		  split_weather = Weather(my_weather.settings,my_weather.simdata,bounds);
%			plot(split_weather,split_weather.t0,split_weather.z0,'wind',1);
%			axis([0 19 0 19]);
%			%axis([0 104500 0 104500]);
%			pause;
%		  new_weather = SubSample(split_weather,time,alt,min_sigma,new_weather);
%
%		  bounds.x0 = p_bar(1)+1;
%		  bounds.xf = my_weather.xf;
%		  bounds.y0 = p_bar(2)+1;
%		  bounds.yf = my_weather.yf;
%		  bounds.z0 = alt;
%		  bounds.zf = alt;
%
%		  split_weather = Weather(my_weather.settings,my_weather.simdata,bounds);
%			plot(split_weather,split_weather.t0,split_weather.z0,'wind',1);
%			axis([0 19 0 19]);
%			%axis([0 104500 0 104500]);
%			pause;
%		  new_weather = SubSample(split_weather,time,alt,min_sigma,new_weather);

		else

			new_sample = my_weather(time,mean(my_weather.x),mean(my_weather.y),alt);
			%wind = [new_sample.u, new_sample.v, new_sample.w];
			%wind = (V(:,2)'+V(:,3)')/2;
			wind = V(:,3)';
			node = dlnode([mean(my_weather.x) mean(my_weather.y)], wind);
			new_weather.append(node);

		end
